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I. INTRODUCTION 

During the last decade, investigation of the detailed structure and properties of small 
sodium clusters attracted a lot of attention (see, e.g., ^ ]3[] and references therein), be- 
cause namely the sodium clusters were used in such important experimental work as the 
discovery of metal cluster electron shell structure [|J and the observation of plasmon reso- 
nances ||, H, [/J. These experiments were definitely among those, which clearly demonstrated 
that atomic clusters and small nanoparticles are in fact new physical objects possessing 
their own properties. The novelty of cluster physics is also greatly connected with the 
fact that cluster properties explain the transition from single atoms or molecules to solid 
state. Comprehensive survey of the field can be found in review papers and books, see, e.g., 
i 1, |, 1, 0,0.0- 

With the discovery of electronic shell structure in free alkali clusters [|], [13|] the essential 



role of the quantized motion of delocalized valence electrons in the mean field created by 
ions in a cluster has been understood. Under different experimental conditions, the detailed 
ionic structure has been found not to affect the properties of alkali and other simple metal 
clusters very much (see, e.g., [§, |], [K| for review). This behavior suggests the validity of 
a jellium model, defined by a Hamiltonian which treats the electrons in the usual quantum 
mechanical way, but approximates the field of the ionic cores by treating them as a uniform 
positively charged background. This naturally leads to a description of the electron density 
in terms of single particle wave functions that extend over the entire cluster. 

Initially, jellium calculations for metal clusters were based on the density functional for- 
malism with the use of pseudopotentials for the description of electron relaxation effects and 
lattice structure [Q. Fully self-consistent calculations for spherical jellium metal clusters 



have been performed within the framework of the spin-density-functional method |L5[ and 
the Kohn-Sham formalism for the self-consistent determination of electron wave functions 
|T3| , [U| . The Hartree-Fock (HF) scheme for the self-consistent determination of the electron 



wave functions of spherical jellium metal clusters was also introduced later in [T7|, [18]. 

Shortly after the discovery of electronic shell structure in free alkali clusters 0, [13] it was 
realized that the detailed size dependence of ionization potentials and other characteristics 
of small metallic clusters can be understood as a consequence of non-spherical cluster shapes 



T9| by analogy with the nuclear shell model. Direct evidence for cluster deformation was 
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achieved in experiments on photoabsorption, where splitting of plasmon resonances caused 
by cluster deformation was observed (see p, |9], |10| and references therein). 

Kohn-Sham calculations for spherical Jl3], [16|] metal clusters have been generalized for 
spheroidal |2(], E_I] and more general axial shapes of light clusters. Light clusters of 
arbitrary shapes have been studied by means of the "ultimate" jellium model p3l Ell. The 



Hartree-Fock approach was generalized for axially deformed cluster systems in [ 25 , EB1 . The 

1 E7|. It was shown that 



existence of the different shape isomers were discussed in |22 



the shape of the magic clusters can deviate from the sphere if higher multipole deformations 
are taken into account. For example, the magic cluster Na^o is not spherical if octupole 
deformations are allowed. This possibility have been pointed out in ]28| and confirmed by 
the 'ultimate' jellium |29j and the Born-Oppenheimer local-spin density molecular dynamics 
method ][3(J. It has been shown that alkali-metal clusters have similar shapes with small 



atomic nuclei |31| . This similarity is a universal result of the density-functional theory. 

The shell-correction method known from nuclear physics |32| has been used in a number 
of papers in studying spheroidal ground state deformations, energetics, stability towards var- 



ious fragmentation channels of metal clusters [33, 34, 35, 36]. The role of higher multipole 



and tri-axial deformations has been elucidated in [3/ 



4ip . Using the shell- 



correction method, the investigation of metal cluster electronic properties such as binding 
energies, ionization potentials, electron affinities, energetics of fission channels and system- 
atic comparison of the theoretical results with the available experimental data have been 
done in |39| (see also j|2| and references therein for a review). 

Dynamical jellium model for metal clusters, which treats simultaneously collective vibra- 
tional modes (volume vibrations, i.e. breathing, and shape vibrations) of the ionic jellium 
background in a cluster, quantized electron motion and interaction between the electronic 



and ionic subsystems was developed in |43], fPJ. This model allowed the widths of electron 
excitations in metal clusters beyond the adiabatic approximation to be described. 

The jellium model provides a very useful basis for studying various collision processes, 
such as photabsorption |45| , photoionization |46|, [47]], elastic |48|, f|9| and inelastic scattering 
jDI O, |52| , electron attachment |53|, I 



| , photon emission pq , pq] , atomic cluster fission 
process WA |58| and others, involving metal clusters. On the basis of the jellium model one 
can develop ab initio many-body theories, such as the random phase approximation with 
exchange or the Dyson equation method and effectively solve many-electron correlation 
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problem even for relatively large cluster systems containing up to 100 atoms or even more. 
Review of these methods in their application to the electron scattering of metal clusters one 



can find in ||59|| . As elucidated in the papers cited above, many-electron correlations are 
quite essential for the correct description of various characteristics of the cluster systems. 

Structural properties of small metal clusters have been widely investigated using quantum 
chemistry ab initio methods. Here we refer to the papers [0, |60], |HT|, ^ |B3|, 0, ^)7| , in 



which optimized geometries, binding energies, ionization potentials, electron structure and 
electron transport properties of small lithium and sodium clusters have been calculated. 

In spite of the fact that both jellium model results and results of ab initio frameworks do 
exist in literature there have been performed no systematic comparison of the results of the 
two different theoretical schemes so far. We fill this gap in our present paper and demonstrate 
that such a comparison is rather illustrative and explains essential physical aspects of the 
formation of various cluster characteristics and properties. Also, we compare the results of 
the Hartree-Fock and local-density approximation (LDA) deformed jellium models and on 
this basis elucidate the role of many-electron correlation effects in the formation of cluster 
deformations. 

On the basis of comparison of the jellium model results with those derived within the 
ab initio theoretical framework and experiment we elucidate the level of applicability of the 
jellium model to the description of alkali cluster properties. For neutral and singly charged 
sodium clusters with N < 20, we have calculated on the basis of Hartree-Fock and LDA 
deformed jellium model the binding energies per atom, ionization potentials, parameters 
of deformation and optimized values of the Wigner-Seitz radii. These characteristics are 
compared with the results derived from the ab initio all-electron theoretical framework for 
the calculation of the cluster ionic and electronic structure based on the density functional 
theory as well as on the post Hartree-Fock perturbation theory on many-electron correlation 
interaction. Comparison performed in our work demonstrates the great role of cluster de- 
formations in the formation cluster properties and the quite reasonable level of applicability 
of the deformed jellium model. 

Our paper is organized as follows. In section II, we provide a brief review of theoretical 
approaches and methods used in the calculation. In section [Til], we present and discuss 
jellium and ab initio results and make their comparison. In section [TV], we draw a conclusion 
to this paper. 
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We use the atomic system of units h = \e\ = m e = 1 in this paper. 



II. THEORETICAL METHODS 



In this work we calculate the binding energies per atom, ionization potentials, deformation 
parameters and the optimized values of the Wigner-Seitz radii for neutral and singly charged 
sodium clusters with N < 20 using the Hartree-Fock and LDA deformed jellium model. The 
jellium model results are compared with those derived in || on the basis of ab initio all- 
electron simulations of the cluster electronic and ionic structure. Below we present a brief 
review of the theoretical methods used in our work. Since the main part of calculations have 
been done within the framework of the deformed jellium model, we focus in our brief review 
on the jellium model approach rather than on the details of the more sophisticated ab initio 
methods for which we refer to ||. 



A. Two-centered jellium model and cluster shape parameterization 

According to the main postulate of the jellium model, the electron motion in a metallic 
cluster takes place in the field of the uniform positive charge distribution of the ionic back- 
ground. Originally, the Hartree-Fock model for metal-cluster electron structure has been 



worked out in the framework of spherically symmetric jellium approximation in [[17], [18[ . 
It is valid for clusters with closed electronic shells that correspond to magic numbers (8, 
20, 34, 40, 58,...). For metal clusters with arbitrary number of valence electrons an open- 



shell two-center jellium Hartree-Fock approximation has been developed (see [25, 26, |57|j ). 
The two-centered jellium method treats the quantized electron motion in the field of the 
spheroidal ionic jellium background, whose principle diameters a and b can be expressed as 
follows: 

R = TgN 1 / 3 is the radius of an undeformed spherical cluster with N atoms, r s is the Wigner- 
Seitz radius, which for the bulk sodium is equal to 4.0. The deformation parameter 5 
characterizes the families of the prolate (5 > 0), and the oblate (5 < 0) spheroids of equal 
volume V c = Aitab 2 /?> = 4irR 3 /3. 
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The electrostatic potential U of the ionic background can be determined from the solution 
of the corresponding Poisson's equation: 



AU(r) = -47rp(r), (2) 



where 



(x 2 + y 2 )lb 2 + z 2 /a 2 < 1 

(3) 

(x 2 + y 2 )/b 2 + z 2 /a 2 > 1 

describes a uniform distribution of the ions charge density in the volume of the cluster. Here 
p c = Z c /V c is the ionic charge density inside the cluster, and Z c is the total charge of the 
ionic core. 

B. Hartree-Fock and LDA formalism 

The Hartree-Fock equations can be written out explicitly in the form (see, e.g., p8| ) 

(-A/2 + U + Uhf) I a >= e a | a > . (4) 

The first term here represents the kinetic energy of electron a, and U its attraction to the 
cluster core. The Hartree-Fock potential Uhf represents the average Coulomb interaction 
of electron a with the other electrons in the cluster, including the non-local exchange inter- 
action, and e a describes the single electron energy. 

According to the density-functional theory, the ground state energy reaches its minimum 
as a function of the density of the system at the exact density ||69|| . A self-consistent method 
for calculation of the electronic states of many-electron systems was proposed by Kohn and 



Sham |[70|| . This method leads to the Kohn-Sham LDA self-consistent equations: 

(-A/2 + U + U H + V xc ) | a > = e a \ a > . (5) 

Here Uh is the Hartree potential, which represents the direct Coulomb interaction of 
electron a with other electrons in the cluster, but does not take into account the non- 
local exchange effects, while V xc is the phenomenological density dependent local exchange- 
correlation potential. The important feature of the LDA method consist in the fact that it 



takes into account many-electron correlations (see, e.g., [71], 72 1 for review). In the present 
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work we use the Gunnarsson and Lundqvist model J73| for the LDA electron exchange- 
correlation energy density e xc , which reads as 

MMr)) = -\ (JL^ - 0.0333 G (r.(r)/11.4) . (6) 

Here r s (r) = (3/47rp e ;(r)) 1/ ' 3 is a local Wigner-Seitz radius, while p e i(r) is the electron 
density in the cluster, and the function G(x) is defined by following relation: 

G(x) = (l + x 3 )ln(l + ^-x 2 + |-i (7) 

The exchange-correlation energy density e xc (p e i(r)) , defines the LDA exchange-correlation 
potential V xc (p d (r)) as 

KeM r)) = S Mpf^m = (8 ) 
9 \ 1/3 1 „ „„„„ , A 11.4 



- 0.0333 In 1 + 



An 2 J r a (r) \ r s (r) 

The Hartree-Fock (||) and LDA (|^) equations have been solved in the system of the 

prolate spheroidal coordinates 74j] as a system of coupled two-dimensional second order 



partial differential equations. The partial differential equations have been discretized on a 
two-dimensional grid and the resulting system of linear equations has been solved numerically 



by the successive overrelaxation method [75] . The third dimension, the azimuthal angle has 
been treated analytically. 

An important characteristic of the cluster, which defines its stability is the total energy 
Etot{N,5). The total energy depends on the size iV of the cluster and its core deformation 
5. The total energy E tot (N,5) is equal to the sum of the electrostatic energy of the ionic 
core E core (N, 5) and the energy of the valence electrons E e i(N, 5): 

E tot (N, 5) = E core (N, 5) + E el (N, 5). (9) 
The electrostatic energy of the cluster ionic core is equal to 

E core (N, S) = \J P(r)t/(r)dr. (10) 
In the HF approximation, the electronic energy E e i(N, S) is given by the general expression 
& 

E% F (N, 5) = < a I ~ A / 2 + U\a> + 

a 

\ I^Qb [c(abk)F k (a, b) + d(abk)G k (a, b)] , (11) 



abk 
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where a and b run over all shells. The values F k (a,b) and G k (a,b) in the Eq. ( |TTD are the 
Coulomb and exchange Slater integrals, q a and are the occupation numbers for orbitals 
a and b, respectively. The Hatree-Fock coefficients c(abk) and d(abk) for the Coulomb and 
exchange energy contributions depend on the occupation numbers (see for details |68|). 
In the framework of LDA the electronic energy of the system is given by [7(| : 

E^ DA {N, 5) = < a I " A / 2 + U\a> + 

a 

1 J P^k^l drdr ' + J p el (r)e xc (Pei(r))dr, (12) 
where the latter term represents the exchange-correlation energy. 

C. All-electron ab initio Hartree-Fock and LDA calculations 

When performing all-electron ab initio Hartree-Fock or LDA calculations for metal clus- 
ters, one has to solve equations which are symbolically the same to the Hartree-Fock (|J) and 
the Kohn-Sham (|5]) equations initially written for the jellium model case. In ab initio calcu- 
lations, instead of using the spheroidal jellium parameterization (|3|) for the distribution of 
the ionic charge density and solving the equations only for the valence electrons, one has to 
deal with the exact Coulomb forces of all the ions and to solve equations (|j) or (|) assuming 
that all electrons are present in the system. The explicit form of the exchange-correlation 
density functional in equations (|5|) can be chosen differently. There are many different func- 
tional forms, although there is no unique one. For example in ||, the ab initio results have 
been obtained with the use of the gradient-corrected Becke-type three-parameter exchange 
functional [|76| paired with the gradient-corrected Lee, Yang and Parr correlation functional 
(B3LYP) [[TT], [77|| . The post Hartree-Fock theories accounting for many-electron correlations, 
such as for example the M0ller-Plesset perturbation theory of the fourth order (MP4) |78j 
have also been used in [|3] for metal cluster simulations. 

In this paper we do not present the explicit forms of these functionals and omit the 
discussion of the methods used for the solution of the Hartree-Fock and LDA equations 
in the ab initio framework. Instead, we refer to the papers, where these functionals and 
methods are presented (see, e.g., || [76], 77, 7J, |7^, 80, 81, 82] and references therein). 
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III. NUMERICAL RESULTS AND DISCUSSION 

In this section we present the results of systematic numerical calculations performed for 
the small sodium clusters in the size range N < 20 on the basis of deformed Hartree-Fock 
and LDA jellium model. We determine the binding energies per atom, ionization potentials 
and cluster deformations. These Hartree-Fock and LDA results are compared with each 
other, with the available experimental data, and with the results derived from the ab initio 
theoretical framework. This comparison elucidates the important role of many-electron 
correlations in metal clusters and establishes the level of applicability of the deformed jellium 
model. 

A. Total cluster energy minimization 

In the Hartree-Fock and LDA jellium models the total cluster energy (|ll]) and ([12]) 
depends on the deformation parameter 5 introduced in (|l|). 

Varying the total cluster energy on 5 one can find its minimum at certain 5 for each 
electron configuration considered. This value of 5 corresponds to the stable geometrical 
configuration of the cluster. 

In figures [1| and we present the total cluster energy per atom E to t(N, 5)/N calculated 
as a function of deformation parameter 5 for the neutral and singly charged sodium clusters 
in the size range N < 20. 

Solid lines in figures [l] and ||] present E tot (N,5)/N for those electronic configurations 
of clusters, that result in the absolute minimum of the total cluster energy. The 5-value 
corresponding to the minimum of E tot (N,5)/N characterizes the cluster shape at the equi- 
librium point. Dashed curves show E tot (N,5)/N for the electronic configurations providing 
the energy minimum the most closely located to the absolute minimum of the total energy. 
These figures demonstrate that for many clusters with open electron shells both oblate and 
prolate isomers are possible and have close energies. The type of deformation that develops 
in the cluster is determined by the type of the corresponding electronic configuration. If 
the electronic orbitals characterizing the chosen electronic configuration are alongated with 
respect to the cluster axis of symmetry then the prolate deformation of the ionic jellium 
background is preferable. In the opposite case, when electronic orbitals lay mostly at the 
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FIG. 1: Total energies per atom of neutral sodium clusters versus deformation parameter 5 calcu- 
lated in the LDA deformed jellium model for different electronic configurations. Solid lines show 
results for those electronic configurations of clusters, that result in the absolute minimum of the 
total cluster energy. Dashed curves show E tot (N,8)/N for the electronic configurations provid- 
ing the energy minimum the most closely located to the absolute minimum of the total energy. 
Cluster images represent correctly the relative cluster sizes, as well as type and value of cluster 
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FIG. 2: The same as Fig. [T], but for the singly charged sodium cluster ions. 
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plane perpendicular to the cluster axis of symmetry the oblate deformation of the cluster 
becomes energetically more favorable. 

It is interesting that the deformed jellium model predicts the existence of a non-spherical 
isomer for the Na>2o magic cluster. So for Nci2o, the second minimum in E to t{N,S)/N 
arises at 5 ~ 0.37, while the electronic configuration keeps the same. In spite of the fact 
that this second local minimum is energetically unfavorable, this prediction is qualitatively 
correct, because it corresponds to the result of ah initio calculations |§, which prove the 
existence of the two cluster isomers (pyramid and deformed structure based on the two 
linked icosahedrons) with the close energies. The pyramid cluster isomer possesses the 
tetrahedral group of symmetry, which is rather high. Thus, this isomer is analogous to the 
spherical cluster in the jellium picture. The deformed icosahedral structure corresponds 
to the prolate jellium cluster configuration. This correspondence would probably be even 
better if one allows the tri-axial deformations in the jellium approach. 

The physical reason for the non-monotonous behaviour of the total energy E to t{N, S)/N 
with increasing the deformation parameter 5 is the strong mixing of the highest occupied 4<r 
state with the unoccupied 5a state. These levels exhibit a rather peculiar behaviour with the 
deformation parameter 5 avoiding each other at 5 ~ 0.37. Such an avoided crossing-point is 
linked to the Wigner's no-crossing theorem ||83|| , which states that two energy levels of the 
same symmetry cannot cross. Thus, the strong interaction of the 4cr state with the more 
prolate type 5<r state results in the appearance of the second minima at 5 ~ 0.37 for Na2o 
cluster. 

Figure |2] demonstrates that similar behaviour of the total energy upon deformation pa- 
rameter observe for singly charged cluster ions. This figure shows that in the jellium ap- 
proach the closed shell cluster ions, Nag and Na^i, have the spherical shape, 5 = 0, at the 
equilibrium point, similarly to the neutral magic clusters Na$ and iVc^o- 

B. Cluster deformations 

In the axially deformed jellium model clusters can either be spherical or have a shape of 
ellipsoid of revolution (spheroid). The spheroidal shape can be of the two types prolate or 
oblate, depending on the sign of the deformation parameter 5 introduced in ([!]). In the ah 
intio approach, the cluster shape is determined by the optimized coordinates of all the ions 
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and it can be characterized by the tensor 

Rij = ^ ^ XjXj (13) 

Here, the summation is performed over all ions in the system. The principle values of this 
tensor R xx , R yy and R zz define the dimensions R x , R y and R z of the ionic charge distribution 
in the cluster along the principle axes x, y and z. Note that tensor R^ is closely connected 
with the cluster moment of inertia tensor and the quadrupole moment tensor of the ionic 
distribution. 

The tensor R^ can also be defined for the jellium model. In this case, sum in (|T3"| ) should 
be replaced by the integral and the integration to be performed over the homogeneous 
spheroidal distribution of the ionic density in the cluster. Then, the principal values of the 
tensor R^ can easily be determined. The result of this calculation reads as 



b 2 a 2 

Rxx = Ryy = jN, R zz = — iV. (14) 

Here, a and b are the principle diameters of the spheroid defined in ([!]). 
In figure [| we present the principle values R xx , R yy and R zz calculated for the neutral 
sodium clusters with N < 20 in the framework of the deformed jellium model according to 



(|14]). The diameters a and b have been determined by minimizing the total cluster energy 
in the LDA approximation (|P2"D as explained in the previous subsection. This calculation 
has been performed with r s = 4.0 (Fig. |^a), which corresponds to the density of the bulk 
sodium, and with optimized value of the Wigner-Seitz radius (Fig. ^b). The cluster energy 
minimization on Wiger-Seitz radius will be discussed in section [111 b\ in more detail. The 
LDA deformed jellium model results are shown in figure [3] by the filled triangles. The filled 
triangles pointing up correspond to R xx =R yy , while those pointing down to R zz . The opened 
triangles are the results of the all-electron ab initio framework derived in |§ with the use 
of the B3LYP density functional. The opened triangles pointing up and down show R xx 
and R zz respectively, while the opened triangles pointing right represent R yy . In figure (| we 
present the results of similar calculations performed for singly charged ions. The notations 
used in figure f| are the same as in figure |3|. 

Comparison of the results presented in figures || and |], demonstrate that the optimization 
of the cluster energy on the Wigner-Seitz radius does not change much the neutral cluster 
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radius optimization (a) and with optimization on r s (b) as a function of cluster size calculated in 
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FIG. 4: The same as Fig. ||, but for singly charged sodium clusters. 
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geometries. For cluster ions, the alteration of the cluster shape with the variation of the 
Wigner-Seitz radius is more noticeable, although it does not improve the agreement of the 
jellium model and ab initio results. This comparison demonstrates that the Wigner-Seitz 
radius variation in the jellium model does not actually improve the quality of the model. 

Figures |3] and f| demonstrate rather good agreement of the jellium model and ab initio 
results. In most of the cases the jellium model predicts correctly the type of the dominant 
cluster deformation, prolate or oblate one. Of course, ab initio calculations include tri-axial 
deformations of the cluster, which turned out to be noticeable for the clusters with the 
open subshells and play important role for clusters like Na\2 — Na±A, Nan, Na^ 3 , Naf 4 , 
iVa^Q. The axially symmetric deformed jellium model does not take into account tri-axial 
deformations and thus in this case always R xx = R yy . 

The axially symmetric jellium model gives the wrong type of deformation in the open 
shell clusters, like Na§, iVa 16 — iVaig, Na§, Naf 7 — Na^Q. However, it is necessary to note 
that for all these clusters there are almost degenerate oblate and prolate isomers within the 
axially symmetric jellium model, as it is shown in figures [I] and ^| and in tables | and |I| (see 
Appendix). Thus, accounting for tri-axial deformations in these clusters plays the crucial 
role as it becomes clear from the comparison of the jellium and ab initio results. 

For the magic clusters Na 8 and Na 2 o, the principle values R xx = R yy = R zz are almost 
identical in both approaches, which demonstrates the closeness to the sphericity of the ab 
initio magic cluster shapes. Note that for the magic cluster ion, Nag, some small defor- 
mation remains in the ab initio approach, while in the jellium approach it turns out to be 
spherical. This demonstrates that the ionic framework of the cluster is not that deformable 
as it follows from the jellium model. 

In Appendix, we compiled in tables the optimized Wigner-Seitz radii r s , total energies per 
atom E tot (N, S)/N, deformation parameters S and the second derivatives of the total energy 
on cluster deformation at the 5-point corresponding to minimum of the total energy. The 
later characteristic of the cluster, d 2 E to t{N,S)/dS 2 , is directly connected to the frequency 
of cluster surface vibrations (see, e.g., for details). All values which are presented in 
the tables [I] and [TT] have been calculated in the LDA deformed jellium model for neutral and 
singly charged clusters. 
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C. Binding energies per atom 

In this paper we calculate the dependence of binding energy per atom in the deformed 
jellium model and compare it with ab initio results from ||. The binding energies per atom 
for the neutral and singly charged clusters are defined as follows: 

E b /N = Ex-En/N (15) 
E+/N = ((N - l)E, + E+ - E+) /N, (16) 

where E^ and E~^ are the total energies of a neutral and singly- charged N-atomic jellium 
cluster respectively. 

Figure [5] shows the dependence of the binding energy per atom for neutral (Fig. |^a) and 
singly charged (Fig. [|b) clusters as a function of cluster size calculated in the deformed 
jellium model. We compare the calculated dependences with the ab initio results from || 
obtained by the B3LYP and MP4 methods. In figure [5] we show the jellium model results 
obtained with bulk and optimized values of the Wigner-Seitz radius. It is seen that the 
cluster optimization on the Wigner-Seitz radius brings the cluster energies down and makes 
them closer to the ab initio results. 

Figure [5] demonstrates that the general trend of the curves calculated within the jellium 
framework turns out to be very close to the one obtained from the ab initio calculation. 
The similarity of the the jellium and ab initio curves is higher for N < 10. In the region 
10 < N < 20 small discrepancy in the behaviour of the curves can be attributed to the 
tri-axial cluster deformations taken into account in the ab initio approach and omitted in 
the axially symmetric jellium model. 

Note that the jellium model results for both neutral and singly charged sodium clusters 
are somewhat closer to the predictions of the MP4 method. This method is based on the 
accounting of the many-electron correlations up to the fourth order of the perturbation 
theory and is free of any adjustable parameters. 

Figure |5] demonstrates that in spite of the simplicity, the jellium model turns out to be 
rather reliable approximation able to reproduce reasonably well the dependence of binding 
energy per atom for both neutral and singly charged sodium clusters. 
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FIG. 5: Binding energy per atom for neutral (a) and singly charged (b) sodium clusters as a 
function of cluster size calculated in the LDA deformed jellium model and compared with ab initio 
B3LYP and MP4 results from §. 
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D. On the role of exchange and many electron correlation interaction. 

To illustrate the importance of the exchange and many-electron correlation interaction in 
the cluster, we plot on figure |6| the total energy per atom calculated both in the HF and LDA 
spherical jellium approximations for neutral (Fig. ^|a) and singly charged (Fig. |6|b) spherical 
sodium clusters as a function of cluster size. The spherical cluster shape was assumed in 
this calculation for the sake of simplicity. Figure |6| shows the significant difference between 
the HF and LDA results for both neutral and singly charged clusters, which is the result of 
the accounting for the many-electron correlation interaction within LDA. 

Figure |^ demonstrates that the LDA total energy dependence possesses the local minima 
at the shell closings N = 2,8, 20 for neutral clusters and N = 3,9, 21 for singly charged ones. 
The HF total energy curve has the extra minima at the half-shell closings, i.e. N = 5, 13, 19 
for neutral clusters and N = 6, 14, 20 for singly charged cluster ions. The appearance of 
these extra minima is the result of the more accurate accounting for the exchange interaction 
within the HF approximation. We found that with increasing cluster size within the given 
shell the favorable electronic configuration of the cluster changes resulting in the formation 
of the extra minima on the HF total energy dependence. Qualitatively, this situation can be 
understood on the basis of the Hund's rule which states that the lowest energy level in the 
system at a fixed electronic configuration is characterized by the maximum value of total 
spin and the maximum possible (at this spin) angular momentum. As an illustration, figure 
|7| shows the ^-dependencies of the total energy per atom for Na^ cluster obtained in the 
Hartree-Fock (Fig. [7(a) and LDA (Fig. ^b) deformed jellium model for different electronic 
configurations. Figure 0a clearly demonstrates that the extra minima on HF total energy 
curve correspond to the electronic states in which spins of all the electrons from the open 
shell are parallel. The LDA framework, with the exchange-correlation potential @, which 
we have used in our work does not take into account the spin polarization effects and thus 
we found no minima on the LDA total energy curve at closing of the half filled shells (see 
Fig. 0b). 
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FIG. 6: Total energy per atom for spherical neutral (a) and and singly charged (b) sodium clusters 
calculated in the Hartree-Fock and LDA spherical jellium model (HF SJ and LDA SJ). 





FIG. 7: Total energy per atom for Na^ cluster versus deformation parameter S calculated in 
Hartree-Fock (a) and LDA (b) deformed jellium model for different electronic configurations. 
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FIG. 8: Ionization potentials of neutral sodium clusters calculated in the deformed jellium model 
and compared with ab initio results from || and with experiment ||. 

E. Ionization potentials 

Another important characteristic of the cluster system is its ionization potential. The 
ionization potential is determined by the energy needed to take an electron out of the cluster. 
It is equal to: 

Vi = E%- E N (17) 

In figure H, we present the ionization potential of neutral sodium clusters calculated 
within the jellium model as a function of cluster size. We compare the jellium model results 
with those obtained in || using ab initio theoretical framework and with the available 
experimental data ||. This comparison demonstrates that the jellium model reproduces 
correctly most of the essential features of the ionization potential dependence on N. Some 
discrepancy, like in the region 11 < N < 14, can be attributed to the neglection of the 
tri-axial deformation in the axially symmetric jellium model. 
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In spite of the fact that ab initio results are closer to the experimental points, one can 
state quite satisfactory agreement of the jellium model results with the experimental data, 
which illustrates correctness of the jellium model assumptions and its applicability to the 
description of sodium clusters. 

Figure |] also demonstrates the role of cluster deformations on the formation of the odd- 
even oscillations in the dependence of the cluster ionization potential on N. Indeed, for 
spherically symmetric clusters this dependence turns out to be monotonous within the range 
of the given shell [f|5|] contrary to the experimental observations. Allowing for the cluster 
deformation and introducing a single deformation parameter 5, we have achieved much better 
agreement of theoretical results with the experimental data as it is clear from figure || 



F. Wigner-Seitz radius variation 

Calculations of the cluster total energy are usually performed at the certain value of the 
Wigner-Seitz radius r s . The bulk value of the Wigner-Seitz radius for sodium is equal to 
4.0. However, one can also perform the calculation minimizing the total cluster energy by 
variation of the Wigner-Seitz radius. 

Figure |9| demonstrates the dependence of the optimized Wigner-Seitz radii on cluster 
size calculated for neutral and singly charged sodium clusters within the HF and LDA 
approximations. This figure shows that the alteration of the optimized r s values is much 
larger for the cluster ions as compared to the neutral clusters. For neutral clusters, the 
optimized values are somewhat larger than the bulk value r s = 4.0 in both LDA and HF 
approximations. The LDA dependence goes closer to the bulk limit. With increasing N 
this dependence approaches the bulk limit, being very close to it also for the magic numbers 
N = 8 and N = 20, which is another manifestation of the shell effect. 

The difference between the optimized values of the Wigner-Seitz radii obtained in the 
HF and LDA approximations can be atributed to the manifestation of the many-electron 
correlation interaction in the system. To illustrate this fact, we plot in figure PD| the difference 



between the optimized Wigner-Seitz radii calculated in the LDA and HF approximations 
for neutral and singly charged sodium clusters. These dependences have the prominent 
peculiarities. The origin of these peculiarities is the same as for those in figure ^ (see section 
HID and discussion therein), although now we consider deformed cluster systems. 
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FIG. 9: Optimized Wigner-Seitz radii for neutral and singly charged sodium clusters calculated as 
a function of cluster size in the HF and LDA deformed jellium models. 

IV. CONCLUSION 

In this paper we performed systematic calculation of various characteristics of neutral 
and singly charged sodium clusters with iV < 20 on the basis of the deformed HF and 
LDA jellium models. We compared the results of our calculations with the ab initio results 
obtained in and with the available experimental data. From these comparisons, we have 
established the level of applicability of the jellium model to the description of various cluster 
characteristics. 

Our consideration shows that the deformed jellium model provides qualitatively correct 
description of the sodium clusters and their ions. The quantitatively reliable results with 
the accuracy below than 10 per cent one can expect from the jellium model description 
providing one allows for the tri-axial cluster deformations. 

We have performed our calculations using HF and LDA approximations. Comparison of 
the results of the two approaches allowed us to demonstrate the importance of the many- 
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FIG. 10: Difference between the optimized values of the Wigner-Seitz radius for neutral and singly- 
charged sodium clusters calculated as a function of cluster size in HF and LDA deformed jellium 
model. 

electron correlations in the formation of cluster characteristics and properties. 

We have performed our calculations for sodium clusters. However, most of the conclu- 
sions should be applicable to other alkali clusters, potassium for example. The level of 
applicability of the jellium approach to other metals, like alkali-earth, requires a separate 
careful consideration. 
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APPENDIX A: TABLES 



In Appendix, we present tables of the optimized Wigner-Seitz radii, total energies per 
atom, deformation parameters and the second derivatives of the total energy on cluster 
deformation calculated in the LDA deformed jellium model for neutral and singly charged 
sodium clusters. 
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TABLE I: Summary of the optimized Wigner-Seitz radii r s , total energies per atom Ef t(N, 5)/N, 



deformation parameters 5 and second derivatives of the total energy on cluster deformation, 
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